Accurate structure factors from pseudopotential methods 



ON 

o 
o 

(N 
Oh' 

o 



J R Trai£| and D M Bird 

Department of Physics, University of Bath, Bath BA2 7 AY, 
(Dated: June, 1999) 



UK 



Highly accurate experimental structure factors of silicon are available in the literature, and these 
provide the ideal test for any ab initio method for the construction of the all-electron charge density. 
In a recent paper [J. R. Trail and D. M. Bird, Phys. Rev. B 60, 7863 (1999)] a method has been 
developed for obtaining an accurate all-electron charge density from a first principles pseudopotential 
calculation by reconstructing the core region of an atom of choice. Here this method is applied to 
bulk silicon, and structure factors are derived and compared with experimental and Full-potential 
Linear Augmented Plane Wave results (FLAPW). We also compare with the result of assuming the 
core region is spherically symmetric, and with the result of constructing a charge density from the 
pseudo-valence density + frozen core electrons. Neither of these approximations provide accurate 
charge densities. The aspherical reconstruction is found to be as accurate as FLAPW results, and 
reproduces the residual error between the FLAPW and experimental results. 

PACS numbers: 78.70.Ck, 71.15.Hx, 71.15.Ap 
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I. INTRODUCTION 

Pseudopotential methods, particularly within the 
framework of total energy plane-wave calculations, are 
extremely powerful for the ab initio description of large 
system of atoms due to their computational efficiency 
and suitability for structural optimisationflj. However, 
they do not yield the correct charge density of the sys- 
tem studied, but a 'pseudo' charge density that does not 
include core electrons and is incorrect close to atomic 
nuclei. This precludes the direct application of these 
methods to the prediction of any properties of the ma- 
terial that depend directly on the charge density, such 
as hyperfine couplings or X-ray structure factors. In this 
paper we reconstruct the all-electron charge density for 
bulk silicon from a pseudopotential calculation using the 
method described by Trail and Bird[2( (hereafter referred 
to as I) , and from this we derive the X-ray structure fac- 
tors. These are then compared with experimental results 
and the results of other theoretical approximations and 
methods to assess the importance of various assumptions 
often made in the calculation of structure factors, and to 
evaluate the success of the reconstruction method. 

Previous methods [H, 0, HI, @] for solving this recon- 
struction problem have relied on the assumption that 
the potential in the core region is spherically symmet- 
ric in order to decouple the differential equations that 
must be solved, and in many cases the charge density 
itself is assumed to be spherical. The method used here 
does not require this to be the case, and we compare the 
structure factors resulting from assuming spherical sym- 
metry to justify the extra effort necessary to develop an 
aspherical reconstruction procedure. The reconstruction 
method itself is based around the embedding approach 
of Inglesfield Q ■ Results from this localised calculation 
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are used to replace the pseudo-charge density where this 
is incorrect, leading to the required structure factors. 

In section |n] a brief summary of the reconstruction ap- 
proach is given - a full description can be found in I. 
We describe how we obtain the structure factors from 
the reconstruction in section ITTT1 and compare the recon- 
struction results with accurate experimental and theo- 
retical structure factors. Rydberg atomic units are used 
throughout the paper. 



II. RECONSTRUCTION METHOD 

The first step in the reconstruction procedure is to ob- 
tain an accurate approximation for the real space single 
particle Green function of the substrate system, in this 
case bulk silicon. We begin with a total energy pseu- 
dopotential calculation performed with a plane-wave ba- 
sis and using the Local Density Approximation (LDA) 
for exchange and correlation. A plane-wave energy cut- 
off of 400 eV is used and 28 Monkhorst-Pack k-points[|] 
are included in the irreducible wedge of the FCC Bril- 
louin zone. These values are more than sufficient to ob- 
tain essentially perfect convergence of the self-consistent 
density and potential, which allows us to attribute any 
errors in our results to the reconstruction procedure. A 
norm-conserving KerkerQ pseudopotential is used, with 
a maximum core radius of 2.0 au. As explained in I 
the method requires r c to be less than half the nearest- 
neighbour atomic separation in the crystal. The resulting 
self-consistent potential is used to obtain a set of eigen- 
states by direct matrix diagonalisation, at 240 k-points 
in the irreducible wedge of the Brillouin zone and with 
a 200 eV plane-wave energy cut-off. Careful tests have 
been carried out to confirm that these values are sufficient 
to provide structure factors with a precision of order ~ 1 
milli-electrons/atom. The spectral representation [lfj is 
used to construct a Green function from the set of plane- 
wave states. This Green function is then used to obtain 
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an embedding potential, which is a term that is added to 
the Kohn-Sham Hamiltonian for the localised core region 
of an atom of interest. The effect of the embedding po- 
tential is to take into account the lattice of pseudo-atoms 
surrounding the chosen atom. The localised embedded 
Hamiltonian is then solved self-consistently (again using 
the LDA) to obtain the Green function of the embedded 
system, from which the charge density in the core region 
can be obtained (see paper I). 

The reconstruction is performed using the same pa- 
rameters as in I, with the embedding radius chosen as 
the 'touching spheres' radius (in this case r s = 2.222 
au). Again, convergence with respect to all parameters 
has been thoroughly checked. The final result we arrive 
at is for the charge density of a single all-electron atom 
embedded in a lattice of pseudo-atoms. This information 
can be used together with the original pseudo-charge den- 
sity to construct an accurate all-electron charge density 
and hence the structure factors for the crystal. 

III. STRUCTURE FACTORS 

Extremely accurate experimental structure factors for 
silicon have been available in the literature for some time 
[Til [H, [H, HU • These results have been used by a number 
of workers to assess the accuracy of parameterised mod- 
els [l5j], FLAPW and other ab initio methods [Hj] and 
Generalised Gradient Approximations to the exchange- 
correlation potential [17|. In view of the accuracy and 
range of data available, both experimental and theoret- 
ical, the reconstructed silicon charge densities are here 
used to construct structure factors for comparison with 
experimental data and the results of FLAPW calcula- 
tions. 



A. Structure Factors from Reconstructed Charge 
Densities 

To obtain the structure factors we require the Fourier 
coefficients of the charge density, 

Ptotal(g) = n I Ptotai(r)e~ lsr d 3 r (1) 
" Jn 

where f2 denotes the unit cell volume and ptotai (r) is the 
total real space charge density, ptotai consists of the orig- 
inal pseudo-charge density between atoms, and the re- 
constructed total charge density within the embedding 
sphere surrounding each atom. Since this integral is a 
linear operation on the charge density it is possible to 
subtract the contribution to the pseudo-density from the 
embedding regions around each atom and add on the con- 
tributions from a reconstruction calculation. This gives 
the expression 

]T[< co " (g)-< elldo (g)] (2) 



where are the position vectors of the atoms in the unit 
cell. The quantities a recon and a pseudo are the Fourier 
integrals of the reconstructed and pseudo-densities re- 
spectively, carried out over the reconstruction sphere sur- 
rounding each atom, and are given by 

«*(g)= / p(r)e-^ r d 3 v (3) 

J \r — Si | <r s 

where r s is the radius of the reconstruction sphere and 
p(r) is the appropriate charge density, reconstructed or 
pseudo. The original pseudo-charge densities are avail- 
able in reciprocal space directly from the plane-wave cal- 
culation, and the reconstructed charge densities are given 
as an expansion in spherical harmonics [2|, which allows 
Eq. © and © to be evaluated. 

For an atom situated at the origin Eq. ([3]) takes the 
form 

a (g) = 4tt V(- i )'y L (g) f 3 pLir^grydr (4) 

L J ° 

where the charge density has been explicitly written as 
an expansion in spherical harmonics, and the identity 

e^ = An^ l Mir)Y£(ci)Y L (r) (5) 

L 

has been used. The radial integral in Eq. (H|) is carried 
out numerically. In our calculations for silicon we choose 
the primitive unit cell, and reconstruct the core region 
of one atom chosen to be at the origin, hence the inte- 
gral in Eq. Q is carried out over a sphere centred on 
this atom. Other atoms within the unit cell must also 
be taken into account, and in the case of silicon there is 
another atom at (|, j, |) related to the origin by an in- 
version symmetry at (g, g, g). The contribution to Eq. 
© from this atom can be derived from the symmetry 
of the unit cell. If the atom at the origin is related to 
an atom at site s by the space group operator {P|s} (de- 
fined by {P|s}/(r) = / (Pr+s)) [HI, where P is a unitary 
transformation, then the integral, a s is 

a s (g) = / e-^ r {P|s}-V(r) ^ (6) 

J|{P|s}- 1 r|<r 3 

By transforming coordinates this reduces to 

a.(g)= / p(p- 1 r)e- Is(r+s) rf 3 r. (7) 

J\r\<r 3 

For silicon the atom at (t, t, t) is related to the atom at 
the origin by the operator {— 1\ (|, t, t)}, an inversion 
followed by a translation. In this case the above expres- 
sion, together with the expansion around the origin in 
spherical harmonics, yields 

a.(g) =47re-**' V(-l)' H)^(g) [ * p L (r)ji(gr)r 2 dr 

L J ° 

(8) 
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where s = j, |). The transformation results in the 
phase factor, and the inversion results in the power of 
(— 1) present in the sum. 

Eq. ((4]) and ([8]) are applied to both the reconstructed 
charge density and the pseudo-density (expanded in 
spherical harmonics) and are then substituted into Eq. 
@ to yield the structure factor as a function of the re- 
ciprocal lattice vector, g. At first it seems a roundabout 
route to calculate the radial expansion of the pseudo- 
density only to convert this back to a reciprocal space 
representation, but this is the most straightforward way 
of replacing the pseudo-charge density with the recon- 
structed charge density in the sphere around each atom. 
One final point is the position of the origin. The coordi- 
nate system used for the reconstruction has the origin on 
one of the silicon atoms in the unit cell (at 43m) , whereas 
the system normally chosen for crystallographic studies 
has the origin at the inversion centre, (3m) [191 ]. Placing 
the origin at the inversion centre gives real structure fac- 
tors, and the origin can easily be shifted to this point by 
introducing an appropriate phase factor into Eq. |T]), or 
simply by taking the magnitude of the complex structure 
factors. 



Comparison with Experimental Results and 
FLAPW Calculations 
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FIG. 1: Difference between static form factors calculated by 
reconstruction and FLAPW methods (circles), and the differ- 
ence between static form factors calculated by pseudo+core 
and FLAPW methods (crosses). 



Before comparison can be made between the theoreti- 
cal and experimental results two further factors must be 
considered. First, the experimentally measured quantity 
(normally given in the literature) is not the Fourier co- 
efficient of the charge density, p(g) , but the form factor 
fhki, which takes into account the lattice structure. This 
is defined as [ItJ 

fhki = Ptotal (g)/ cos ((h + k + I) ~) (9) 

where (hkl) are the indices of the reciprocal lattice vector. 
For (hkl) values that satisfy the criteria h + k + l = 4n + 2 
for n integer, the denominator on the right hand side is 
zero, and the structure factor is given. 

The second effect that must be taken into account 
when correlating the theoretical and experimental results 
is the thermal motion of the lattice. The majority of ex- 
perimental data for structure factors are taken at room 
temperature, and the thermal energy 'smears out' the 
charge density, reducing the amplitude of the higher or- 
der structure factors. This can be described by a con- 
volution integral in real space, which corresponds to a 
further correction factor in reciprocal space to give the 
dynamic structure factor 

fH y M=hkie- Bg2 ^ 2 (10) 

where B is the Debye- Waller parameter [IH, [l7| • 

Structure factors obtained from the core reconstruc- 
tion are here compared with those obtained from three 



sources: from the simple addition of free atom core states 
to the original pseudo-charge density [l6| , structure fac- 
tors obtained using the FLAPW method by Zuo et al 
[ItJ and structure factors determined experimentally by 
Cumming and Ha rt [ill and Saka and Kato [TJ], as 
given by Zuo et al [IT]]. The pseudo+core structure fac- 
tors are obtained from the charge density of the original 
pseudopotential calculation together with the core charge 
densities of the original atomic calculations used to create 
the pseudopotential. The contribution from the atomic 
core charge density is included at the atomic sites in the 
same manner as described above, ie 

Ptotalis) = Ppseudo(g) + i a sr e (§) ( n ) 

i 

where a c ° re is the contribution from the core states at site 
This structure factor is expected to show significant 
error, since the valence charge density will be incorrect 
close to the atomic sites. 

Zuo et al [TtJ calculated Si structure factors using the 
FLAPW method, and they give results using the LDA, 
and two different GGAs. Since the core reconstruction 
calculation carried out here employs the LDA, the re- 
construction results are only compared with the LDA 
FLAPW results. For a successful reconstruction scheme 
we would expect to accurately reproduce these results, 
since the same physical approximations have been made 
even though the algorithmic implementation of the two 
methods is entirely different. 

We begin by comparing the theoretical form factors, 
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uncorrected for temperature, and only for the (hkl) val- 
ues for which experimental data are available (experi- 
mental data are given in Table HJ. Fig. Q] shows the 
difference between the reconstructed and FLAPW form 
factors, and the difference between the pseudo+core and 
FLAPW form factors, ft is apparent that the reconstruc- 
tion agrees very well with the FLAPW results - the av- 
erage absolute difference for the reconstructed results is 
only 3 milli-electrons/atom whereas for the pseudo+core 
result the average absolute difference is over 25 times 
greater at 76 milli-electrons/atom. 

C. Experimental, FLAPW and Reconstructed 
Structure Factors 

In order to compare the static structure factors given 
above with experimental data a value for the Debye- 
Waller parameter in Eq. (|I0|) is required. This is com- 
monly taken to be a free parameter and varied to min- 
imise the error between the experimental and theoretical 
results. The value of B used here is that employed by Zuo 
et al[I3, B = 0.4668 A 2 . In their paper values of B are 
obtained by minimising the error of high |g| values only, 
for a number of different ab initio methods. These high 
|g| structure factors depend largely on the core states of 
the atoms that make up the lattice, so the best values 
should result from methods that accurately describe the 
core states. Zuo et al found that a calculation of these 
high order structure factors usin g th e Multi Configura- 
tion Dirac-Fock (MCDF) method [2Cj gives the best fit at 
high |g|, and therefore took the associated B parameter 
to be the best estimate. It should be noted that a bet- 
ter fit between experiment and theory can be obtained 
for a different B value, but this would effectively use the 
description of a physical effect, the thermal smearing, to 
adjust for deficiencies in the theory, such as the LDA. 

Table U gives the experimental data, with recon- 
structed, FLAPW and pseudo+core dynamic form fac- 
tors. The quality of the theoretical data is assessed by 
two statistics - the i?-factor and GOF parameter. The 
i?-factor is given by 

Ei rtheory rexpi 
i \ J i J i f-i n\ 

Eii/n [ ' 

and the goodness of fit parameter by 
1 N 

GOF = ± ^l^)Ul he ° rv - ID 2 (13) 

i=l 

where of is the sample variance of the i form factor. 
The variance of is taken to be the average of the esti- 
mated error for all data points in line with the approach 
of Zuo et al, and takes the value (0.0022) 2 /e 2 atom~ 2 . 

From the data in Table U it can be seen that the re- 
construction calculation describes the experimental data 
as well as the FLAPW results. For both sets of data 
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TABLE I: Dynamic form factors from experiment, reconstruc- 
tion, FLAPW and Pseudo+core calculations. Estimated er- 
rors of the experimental data are given in parentheses, and 
the experimental and FLAPW data are taken from Zuo et al 
[l7l ]. A Debye- Waller parameter of B = 0.4668 A 2 is used, 
as calculated by Zuo et al [l7j . 

the i?-factor is 0.24 %, and the GOF is - 35, with the 
GOF for the reconstruction slightly greater than that 
for FLAPW. The average absolute error \f theor v - f ex P\ 
is 10 milli-electrons/atom for both the FLAPW and re- 
construction calculations and 70 milli-electrons / atom for 
the pseudo+core results. The maximum error is roughly 
~ 20 milli-electrons/atom for the FLAPW and recon- 
struction results, and ~ 100 milli-electrons/atom for the 
pseudo+core results. Fig. [5^ shows the residual error 
(Sf = f theor y - f ex P) of the FLAPW results together 
with the error bars of the experimental data, and Fig. 
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FIG. 2: Residual error 

^theory _ jexp^ q{ ^ FLAPW and ( b ) reconstructed dynamic form factors, with error bars of experi- 
mental data shown. Note the change of scale from Fig. [1] 



[2j} the residual error for the reconstruction calculation. 
The errors are very similar, even to the point of a sig- 
nificant correlation existing between the two. This sug- 
gests that the errors present are largely due to the the- 
ory shared by the calculations, specifically the LDA. It 
should also be noted that the data presented by Zuo et 
al is calculated for a lattice constant of ao = 5.4307 A, 
whereas the reconstruction calculations are carried out 
for a a = 5.4300 A. 

Finally we give the i?-factor and GOF parameter com- 
paring the pseudo+core and reconstructed results with 
the FLAPW results. The pseudo+core form factors give 
a i?-factor and GOF of 1.55 % and 1349 respectively, 
while the reconstruction gives 0.06 % and 3.6. 

D. Spherical Symmetry 

One of the strengths of our reconstruction method is 
that it does not require spherical symmetry of the charge 
density in the reconstruction region near the cores of 
the atom. To assess the importance of the aspherical 
components of the charge density we replace the origi- 
nal pseudo-density in the reconstruction sphere with the 
spherical part only of the reconstructed density. Fig. [3] 
gives the residual error of the reconstructed form factors 
from the experimental data for this case. The i?-factor 
is 0.64 % with a GOF of 485 - considerably worse than 
for the full aspherical reconstruction. From this data it 
is apparent that the aspherical components of the charge 
density are essential for the calculation of accurate form 



factors. However, it is interesting to note that if we re- 
place the spherical part of the pseudo-density with the 
spherical part of the reconstructed density (that is the 
charge density components pl for I > within the em- 
bedding sphere are given by the pseudo-charge density) 
we obtain form factors that are almost as accurate as the 
FLAPW and fully aspherical results. In this case the R- 
factor is 0.25 %, the GOF is 37 and the mean absolute 
error is ~ 11 milli-electrons/atom). 

IV. CONCLUSIONS 

In this paper all-electron states have been recon- 
structed successfully from a total energy pseudopoten- 
tial calculation, giving an accurate charge density in the 
region near atomic sites. This reconstruction is carried 
out using the embedding method described in a recent 
paper [2j. The reconstruction calculation itself uses a 
scalar relativistic description for the valence states, in a 
fully aspherical potential and using LAPW basis func- 
tions. The core states are calculated fully relativistically 
by direct solution of the Dirac equation in a spherical 
average of the self consistent potential. It is apparent 
that the reconstruction method itself has a lot in com- 
mon with FLAPW methods. 

Structure factors have been derived from the recon- 
structed silicon charge density for comparison with accu- 
rate experimental data and FLAPW calculations. Agree- 
ment is excellent with both the FLAPW and recon- 
structed form factors agreeing with experimental results 
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with an average absolute error of 10 milli-electrons/atom 
while the experimental data itself is accurate to 3 — 5 
milli-electrons/atom. The FLAPW and reconstructed 
form factors agree extremely well with each other, with 
an average absolute difference of 3 milli-electrons/atom. 
In addition to this the residual errors for both methods 
of calculation show significant correlation, indicating that 
they arise from the physical approximations common to 
both methods. 
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